mm 


NASA Contractor Report 178179 
ICASE REPORT NO. 86-60 

ICASE 


GLOBAL COLLOCATION METHODS FOR APPROXIMATION 
AND THE SOLUTION OF PARTIAL DIFFERENTIAL EQUATIONS 



FCB APPROXIMATION AND TEE SOLUTION OF 
PARTIAL DIFFERENTIAL EQUATIONS Final Report 
(NASA) 68 p CSCL 12A Unclas 

G3/64 44364 


Contract Nos. NAS1-17070, NAS1-18107 


September 1986 


INSTITUTE FOR COMPUTER APPLICATIONS IN SCIENCE AND ENGINEERING 
NASA Langley Research Center, Hampton, Virginia 23665 

Operated by the Universities Space Research Association 


NASA 

National Aeronautics and 
Space Administration 

Langley Rese ar ch Center 

Hampton, Virginia 23665 





GLOBAL COLLOCATION METHODS FOR APPROXIMATION 
AND THE SOLUTION OF PARTIAL DIFFERENTIAL EQUATIONS 


A. Soloraonof f 
ICOMP 

NASA Lewis Research Center 


E. Turkel 

Institute for Computer Applications in Science and Engineering, 

NASA Lewis Research Center 
and 

Tel-Aviv University 


ABSTRACT 

We apply polynomial interpolation methods both to the approximation of 
functions and to the numerical solutions of hyperbolic and elliptic partial 
differential equations. We construct the derivative matrix for a general 
sequence of the collocation points. The approximate derivative is then found 
by a matrix times vector multiply. We explore the effects of several factors 
on the performance of these methods including the effect of different 
collocation points. We also study the resolution of the schemes for both 
smooth functions and functions with steep gradients or discontinuities in some 
derivative. We investigate the accuracy when the gradients occur both near 
the center of the region and in the vicinity of the boundary. The importance 
of the aliasing limit on the resolution of the approximation is investigated 
in detail. We also examine the effect of boundary treatment on the stability 
and accuracy of the scheme. 
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1. INTRODUCTION 

In this study we consider the accuracy of the pseudospectral 
approximation both for a function and also for the numerical solution of 
differential equations. We shall only consider collocation methods, but most 
of the results shown also apply to Galerkin methods. We approximate the 
function, f(x), by a polynomial, p^(x), that interpolates f(x) at N + 1 
distinct points xo,...,x^. f'(x) is approximated by Pjj'(x) which is 
calculated analytically. In solving differential equations we use an approach 
similar to finite differences. Thus, all derivatives that appear are replaced 
by their pseudospectral approximation. The resultant system is solved in 
space or advanced in time for time dependent equations. Hence, for an 
explicit scheme, nonlinearities do not create any special difficulties. 

This approach is equivalent to expanding f(x) in a finite series of 
polynomials related to Xq,...,x^. For a Galerkin method, the coefficients of 
this series are obtained from the infinite expansion. For a collocation 
method, the coefficients are obtained by demanding that the approximation 
interpolate the function at the collocation points. This requires 0(N Z ) 
operations. For special sequencies of collocation points, e.g.,^ Chebyshev 
methods, this can be accomplished by using FFT^s and only requires O(NlogN) 
operations. Every collocation method has two interpretations: one in terms 
of the collocation points and one in terms of a series expansion. In the 
past, this has lead to some confusion. As an example we consider the case of 
a Chebyshev collocation method with x^ = cos(ttj/N). From an approximation 
viewpoint, we know [11, 15 - 18] that the maximum error for interpolation at 
the zeroes of T N (x) is within (4 + 2/^logN) of the minimax error and 
converges for all functions in C^. The bound for the error based on the 
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points x j , given above, is even smaller than this [13]. There also exist 
sharp estimates in Sobolev spaces [3]. Since the minimax approximation has an 
error which is equi-oscillatory we expect the Chebyshev interpolant to be 
nearly equioscillatory. Indeed, Remez suggests using these Xj as a first 
guess in finding the zeros of f - P N in his algorithm for finding the 

minimax approximation. Thus, we would expect that when used to solve 
differential equations that the error would be essentially uniform throughout 
the domain. 

On the other hand, viewed as a finite difference type scheme, one expects 

the scheme to be more accurate near the boundaries where the collocation 

points are clustered. At the center of the domain the distance between points 

is approximately tt/ 2N while near the boundary the smallest distance 

2 2 

between two points is approximately tt /2N . Hence, the spacing at the 
center is about tt / 2 times coarser than an equivalent equally spaced 
mesh. Near the boundary the Chebyshev points are about 4 n/tt times finer 
than an equally spaced mesh. From this point of view, we expect the accuracy 
and resolving power of the scheme to be better near the boundaries. However, 
the bunching of points near the boundaries only serves to counter the tendency 
of polynomials to oscillate with large amplitude near the boundary. We shall 
also consider collocation based on uniformly spaced points. Since, we 

consider polynomial interpolation on the interval [-1,1] we get qualitatively 
different results than obtained by Fourier or finite difference methods even 
for the same collocation points. In fact, we shall see that the boundaries 
exert a strong influence for this case similar to the interpolation based on 
Chebyshev nodes. 

Connected with this, we shall examine the influence of boundary 
conditions on the accuracy and stability of pseudospectral methods. In 
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general , global methods are more sensitive to the boundary treatment than 
local methods. We also consider the effect of the location of the collocation 
points on both the accuracy and stability of the scheme and its effect on the 


allowable time step for an explicit time integration algorithm. 


2. APPROXIMATION AND DIFFERENTIATION 

We assume that we are given N + 1 distinct points x Q < xj < ... < x N . 
Given a function f(x) it is well known how to approximate f(x) by a 

polynomial P N (x) such that P N (xj) - f(xj), j = (),..., N. We define a 

function e k (x) which is a polynomial of degree N and e k (xj) = 
Explicitly, 


e k (x) = T n < x " x » ) 
k 1-0 % 


( 2 . 1 ) 


a k = 11 (x k “ x » >• 
£=0 % 


(2.1b) 


Then the approximating polynomial is given by, 


P N (x) = l f( x k )e k (x), 
N k=0 


( 2 . 2 ) 


We next consider an approximation to the derivative of f(x). We 
construct this approximation by analytically differentiating (2.2). The value 
of the approximate derivative at the collocation points is a linear functional 
of the value of the function itself at the collocation points. Hence, given 
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the N + 1 values f(xj) we can find the values P N (x) 


by a matrix 

multiplying the original vector (f (xq) , . . . ,f (xjp) . We denote this matrix 

* 

by D = (d.,). By construction, DP = P is exact for all polynomials of 

J K. N 

degree N or less. In fact, an alternative way of characterizing D is by 
demanding that it give the analytic derivative for all such polynomials at 
the N + 1 collocation points. In particular, we shall explicitly 

construct D by demanding that, 


De k ( x j) = e k< x j)» j ,k = 


(2.3) 


i.e. , 



- k-th row = 


w 


\ ( V J 


Performing the matrix multiply, it is obvious that 


de k 

d., = -r-=. (x.) 

jk dx j 


(2.4) 


We next explicitly evaluate d ^ in terms of the collocation points 
x j . Taking the logarithm of (2.1) we have 
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N 


log(e ) = l log(x - x ) - log(a ). 


-\r ' L O'" "} • 

1=0 
J l*k 


Differentiating, we have 


N 

e (x) - e, (x) l l/(x - x ). (2.5) 

k £=0 
Jl^k 

In order to evaluate (2.5) at x * Xj , j £ k, we need to eliminate the zero 
divided by zero expression. We therefore, rewrite (2.5) as 


N 


e fc (x) = e fc (x)/(x - Xj) + e fc (x) £ l/(x - x fc ). 

A “0 
1*3 ,k 


Since, e k (xj) =0 for j * k we 


have that 


e k (x.) = lim e k (x)/(x - x. ). 

^ x+x . J 

J 


Using the definition of e k (x), (2.1), we have 


d. t -JL 

jk a 


N 

n 

k 1=0 


a . 
3 


- V - 


V 7 


(see (2.1b)) 


(2.7) 


Af j » k 

While the above formulas (2.5), (2.7) are computationally useable, it is 
sometimes preferable to express the formulas slightly differently. We, 
therefore, rederive these formulas using a slightly different notation. 
Define, 


N 

^ m. i ( x ) = n (x - x ). 
N 1 1=0 


( 2 . 8 ) 
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Then, 


and so 


Vi <x) 


N 


I 

k=0 


N 

n 

£=0 

Ifk 


(x - x t ) 


(2.9) 


4* 


N+l 



(2.10) 


It can also be verified that 


Hence , 



N+l 


( V 


= 2a. 


N 

l ~ 

k=0 X j 

kti 


- X. 


V x j • v 


N i 

l Z 

1=0 x k X Jt 
A=k 


♦ 

*£i 


* 

N+l^ X j ^ 
( x k )(x.- x k ) 


2 ^+i (x k> 


( 2 . 11 ) 


j * k 

(2.12a) 
j = k (2.12b) 


Given a^ it requires another 2N^ operations to find the off diagonal 
elements by (2.12a). It requires operations to find all the diagonal 

elements from (2.12b). Hence, it requires about 4N^ operations to construct 
the matrix D. We multiply the matrix D on the left by 

diag( 1/a^ , . . . ,1/a^) and on the right by the matrix diag(a^ , . . . ,a N ) . 
Then D is similar to a matrix Dj where is a sum of an antisymmetric 

matrix and a diagonal matrix. Since is a polynomial of degree N - 1 
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^N+l cannot be zero at all the collocation points. Hence, the diagonal 
portion of is nonzero. 


In many 

cases 

x 0 = -1, 

x N = 1 and 

the other 

x i 

are 

zeros of 

some 

polynomial 

Q N -i( x ) 

. Hence, 

W x) ■ 

(x 2 - 1)Q n . 

^(x) 

• 

One can 

then 

rewrite the 

formula 

for d Jk 

in terms of 

Qn-i( x ) • 

For 

j. 

k t 0, N 

we 


reproduce the formulas of Tal-Ezer [20]. He further points out that if 
Qn_i( X ) is a Jacobi polynomial associated with the weight function 
(1 - x) a (1 + x) S then 


v 1 _ -(a + 1) _ (g - 1) 

I x k - x^ Ax k - i) 2ix k + 1) 


(2.13) 


where the sum is over the roots of Q n _j(x). This can then be used to 

simplify (2.12b). When the ends points x = -1 or x = 1 are included in 

the collocation points then these must be explicitly accounted for to find 


d jk* 

For the standard Chebyshev collocation points, we have 


Xj = cos(Trj/N) j = 0,...,N. 


(2.14) 


Note that this orders the points in reverse order from our usual assumption. 
In this case one can evaluate the derivative by using FFT's. This requires 
only NlogN operations rather than the operations required by a matrix 
multiply. Computationally it is found that for N < 100 that the matrix 
multiply is faster than the FFT approach, see, e.g., [23]. The exact 
crossover point depends on the computer and the efficiency of the software for 
computing FFTs and matrix multiplies. The matrix multiply has the advantage 
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that it is more flexible and vectorizable. Thus, for example, both the 
location and the number of the collocation points is arbitrary. In order to 
use the FFT approach, it is required that the collocation points be related to 
the Fourier collocation points, e.g. , Chebyshev. Furthermore, the total 
number of collocation points needs to be factorizable into powers of 2 and 3 
for efficiency. The efficiency of these factors depends on the memory 
allocation scheme of the computer. Other collocation nodes than (2.14) are 
considered in [3, 7, 13] • The matrix D for the Chebyshev points (2.14) is 
given in [7 ] . 

In Appendix A, we consider the problem when we have N collocation nodes 
but wish the derivative matrix to be exact (in least squares sense) for M > N 
functions which need not be polynomials. 


3. PARTIAL DIFFERENTIAL EQUATION 

We consider in this study three applications of collocation methods: (1) 
approximation theory, (2) hyperbolic equations, and (3) elliptic equations. 
For approximation theory we need only discuss accuracy. We first need some 
way to measure the approximation error that can be used on a computer. We 
cannot use the error at the collocation points since, by construction, this 
error is zero. Instead, we use 


" f Pn " jJq P N^ X Jl^ 


(3.1) 


for some sequence of points 
general, we shall choose N' 


Xj which are not the collocation points. In 
much larger than N, the number of collocation 
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points. If the original points are chosen as Chebyshev nodes, then we again 
choose the x j , as Chebyshev nodes based on this larger number, N'. Because 
of this selection of nodes the sum in (3.1) approximates the Chebyshev 
integral norm, i.e., 

0 1 [f(x) - P (x)] 2 

dx (3.2) 

17^7 

! 1 £ - 0,N 

2 £ = 0 ,N 

When the collocation points are evenly spaced then we shall choose the nodes 
of the integration formula to be also uniformly spaced. In this case 


N - 


J 

-1 


where 


1 TT 
w rt = — “ 

* C A N 


~ and 


nf - P N II 2 ~ f [f(x) - P N (x)] 2 dx 


where 


w„ 


1 1 
c, ' » • 


(3.3) 


For general collocation points it is not clear how to choose the weights 
w^ in the norm. An alternative possibility is to measure the error in some 
Sobolev norm. In this case, the finite sum can be based on the original 
collocation points and the norm is the L 4 norm of the derivative. In this 
study all errors will be given by (3.2) regardless of the distribution of the 
collocation nodes. 

For hyperbolic problems we need to be concerned with stability in 
addition to accuracy. We also study the influence of the boundary treatment 
on both the accuracy and stability of the method. For simplicity we shall 
only consider the model equation 
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ll = a(x) u x -1 < X < 1, t > 0. (3.4) 

If a(x) is positive at both boundaries then we need to impose a boundary 
condition at x = 1. If a(-l) is negative while a(l) is positive, then we 
impose boundary conditions at both ends. On the other hand if a(-l) is 
positive while a(l) is negative then no boundary conditions need be given. 
For spectral methods, it is important that this distinction be preserved at 
the approximation level. Thus, whenever analytic boundary conditions are not 
given the spectral technique will be used to advance the solution at the 
boundary. The given boundary conditions are always chosen so that we know the 
analytic solution. 

We will solve the differential equation (3.4) by a pseudo-spectral 
algorithm. Thus, we will consider the solution only at the collocation 
points. We then replace the derivative in (3.4) by a matrix multiply as 
described in section 2. We next multiply a(x) at each collocation point by 
the approximate derivative at that point. We now have a system of ordinary 
differential equations in time. To advance the solution in time we could use 
any ODE solver. In particular, we shall use a standard four stage fourth 
order Runge-Kutta formula. This formula has several advantages. First, since 
it is fourth order in time (for both linear and nonlinear problems), it is 
closer to the high spatial accuracy of the spectral method than a second order 
formula. Also, the region of stability includes a significant portion of the 
negative real half plane and so is appropriate for Chebyshev methods which 
have eigenvalues in the negative half plane. Finally, if we look along the 
imaginary axis it has a comparatively large stability region. An alternative 
method is to use a spectral method in time. However, it is difficult to 
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generalize such methods to nonlinear problems while Runge-Kutta methods extend 
trivially to nonlinear problems. 

Since this is an explicit method (even though all the points are 
connected every time step) it is easy to impose boundary conditions after any 
stage of the algorithm. Whenever we wish we can let the pseudospectral method 
advance the solution at the boundary also. Since the method is explicit, 
there is a limit on the allowable At because of stability 

considerations. Heuristically, one can consider this stability limit as 
arising from two different considerations. One is based on the minimum 
spacing between mesh points, which usually occurs near the boundary. As noted 
above, this is heuristic since the domain of influence of each point is the 
entire interval. Alternatively, one can derive a stability limit by finding 
the spectral radius of a(x) times the derivative matrix. This is also 
heuristic since the derivative matrix is not a normal matrix. For a(x) 

constant both methods indicate that At varies with 1/N . The exact 

constant varies with the particular Runge-Kutta method used. For a two stage 
Runge-Kutta method, the stability limit is about three times the minimum 
spacing. For further details, the reader is referred to [5, 7] and the result 
section. 

In Appendix B, we present the proof of the stability of Chebyshev 
collocation at points (2.14) for u t = u x . 

For our model elliptic problem, we shall choose a Poisson equation 


Au = f (x,y) , -1 < x,y < 1 


(3.5) 


- 12 - 


with u(x,y) prescribed on all four sides. As before f(x,y) will be chosen 
so that we know the analytic solution. 

We solve a time independent equation since it is easier to distinguish 
the resolving power of the scheme in different regions of the domain. For a 
time dependent equation, it can be difficult to distinguish local accuracy 
since inaccuracies propagate from one part of the domain to another. This is 
especially true for systems of hyperbolic equations with characteristics 
travelling in each direction. When the time independent equations is 
elliptic, then the solution is smooth. In particular, u(x,y) has at least 
two derivatives even if f(x,y) is only continuous. The smoother f(x,y) is 
the smoother u(x,y) will be, assuming the boundary conditions are 
sufficiently smooth. 


4. RESULTS 

In this section, we describe the computational results that illustrate 
many of the properties of pseudospectral methods. We begin with the 
approximaton of functions. Unless otherwise noted, the collocation points 
will be the Chebyshev nodes, (2.14). As is well known, interpolation at these 
points yields a maximum error which is not much worse (O(logN)) than the best 
possible minimax approximation [13 - 18]. Nevertheless, we shall see that the 
quality of the approximation can vary greatly for different functions. We 
shall also see the effect of varying the collocation points. 

In Figure la, we display the pointwise error in approximating the 
function u(x) = sin(20x-m) where m varies between 0 and tt/ 2. Thus, 
u(x) varies between a sine and a cosine function. The top graph in Figure 1 
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is the error for an approximation to a sine wave. The phase changes in the 
following graphs and the bottom graph is the error for a cosine function. In 
this case we chose 28 Chebyshev collocation points. For m = 0, i.e., a sine 
function, the amplitude of the error is larger. This occurs since sin(x) is 
an odd function and hence the coefficient of T N is zero and so in essence we 
are only using 27 polynomials. This is verified in Figure lb by using N = 
29; for this case the error of the cosine function is larger. Nevertheless, 
this result is interesting for time dependent problem where the solution 
varies between a sine and cosine function. In addition, we also notice that 
for m = 0 the largest errors occur in the middle of the domain while for 
m = tt/ 2 the larger errors are near the boundaries. Thus, for smooth 

functions the maximum error can occur anywhere in the domain. There is no 
need for the error to be smaller near the boundaries where the collocation 
points are bunched together. 

In Figure 2 we show the pointwise error in approximating the function 

u(x) = |x - xq| which has a discontinuous derivative at x = x . We define 

0 

a point as being half way between two nodes in the Chebyshev sense when 

x = cos (ir( j + V 2 )/N) . 

j + V 2 

The top of the graph displays the error when the discontinuous derivative is 

located halfway between nodes while the center of the graph shows the error 

when the discontinuity in the derivative occurs at a node. The other graphs 
show progressively other locations of x j . Thus, we see that when the 

discontinuous derivative occurs half-way between nodes in the Chebyshev sense, 
then the error has a sharp peak near the discontinuity but is close to zero 
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elsewhere. When the discontinuity occurs near a node then the error is more 
spread out and several peaks may occur but the maximum error is decreased. 
Gottlieb has observed similar phenomena in other problems. For other values 
of x the error goes smoothly between these extremes. 

In Figure 3a, we examine the effect of the aliasing error in the 

approximation of a function. Gottlieb and Orszag [5] show that one needs at 

least tt points per wave length when using a Galerkin Chebyshev 

approximation. In Figure 3, we approximate sin(Mrrx) with N Chebyshev 

nodes in a pseudospectral approximation. We plot the error, (3.2), as a 

function of the number of points past the aliasing limit. As before the error 

begins to decrease exponentially when there are tt points per wave 

length. We further see that in order to reach a fixed error the number of 

collocation points, N, should vary (approximately) as the aliasing limit 
1/3 

plus M A/ . Computationally, it is hard to find the exact exponent, but it 

seems to be between 0.3 and 1/3. In Figure 3b, we see that for f(x) = 

tanh(mx) there is no sudden aliasing limit. Rather there is a gradual 
reduction in the error as N increases. 

For a Fourier method, it can be shown that one only needs two points per 
wave length rather than tt points per wave length. It might be speculated 

that this is due to the larger spacing of the Chebyshev method near the middle 

of the domain. In fact, asymptotically, the largest spacing between Chebyshev 
node is exactly tt/2 times as large as for Fourier nodes. In Figure 4, we 
consider the same case as in Figure 3, but where the collocation points are 
evenly spaced. One sees that one again need about it points per wavelength 
before exponential accuracy occurs even though the spacing is the same as for 
the Fourier method. There is a theorem that interpolation based on uniformly 
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spaced points converges for analytic functions. Nevertheless, we see in 
Figure 4 that the approximation begins to diverge if N is sufficiently large 
with respect to M. The calculations for these case were carried out on a 
CRAY computer which has about 15 significant figures. Using double precision 
(about 30 significant digits) one stabilizes the procedure until larger N 
are reached at which point the approximation again diverges. Hence, even 
though the function is analytic nevertheless roundoff errors eventually 
contaminate the approximation. Hence, collocation based on uniformly spaced 
node is risky even for analytic functions because of the great sensitivity of 
these collocation methods to any noise level. 

In Figure 5a, we study the resolving power of Chebyshev methods when 
there are sharp gradients. It is often stated, that Chebyshev methods are 

ideal for boundary layer flows since they naturally bunch points in the 

2 

boundary layer. In Figure 5a, we plot the L error when we are 

approximating the function u(x) = tanh(M(x - Xq)), for M = 8, 32, 128, 512, 
and 2048 and N = 31. As M increases the gradient becomes steeper and in 
the limit approaches a Heaviside function. Furthermore, tanh(x) is also a 
solution to Burger's equation and so appropriate to model boundary layers. We 
see that, indeed, for moderate values of M the accuracy is greater when the 
gradient occurs closer to the boundary. Thus, given a moderate slope a 

Chebyshev collocation method "sees 11 the gradient better if it is near the edge 
of the domain. Thus it may be advantageous to consider multidomain approaches 
[9, 10]. However, when the slope becomes too large so that it is not resolved 
by the collocation points, then the error is equally large everywhere. In 

particular, a true discontinuity, e.g., a shock, is not resolved any better 
near the boundary then it is in the middle of the domain. 
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In the previous case it was implied that the Chebyshev method resolves 
gradients better near the boundary because the nodes are closer together near 
the edge. To check this hypothesis, we plot the same case in Figure 5b where 
now the collocation is based on uniformly spaced points. We consider the same 
case in Figure 5a but now choose M = 2 , 4 , 8 , 16 , 32 . We choose lower values 
of M then before since the approximation based on evenly spaced points does 
not converge when the gradient is too large. For the same M the errors are 
much larger for the uniformly spaced nodes than Chebyshev spaced nodes. 
Nevertheless, the errors are much smaller when the gradients occur near the 
boundary. Thus, gradients in the "boundary layer" are better resolved than in 
the center of the domain even though we are using interpolation based on 

o 

uniformly spaced points. In fact, the ratio of the Ir error when the 

p 

gradient is at the center to the L error when the gradient is near the edge 
is even larger for uniformly spaced nodes than for a Chebyshev distribution of 
nodes. In both cases, we used the Chebyshev norm (3.2). However, the results 
do not depend on the details of the norm. 

In order to explain this phenomenon we examine the singularity of the 
function in the complex plane. To simplify the discussion we shall consider 
the easier case of an expansion of a function in Chebyshev polynomials. In 
this case it is known [14] that the approximation converges in the largest 
ellipse with foci at +1 and -1 that does not contain any singularities. 
The equation of an ellipse with foci at ±1 is 



= 4 


where 


> 2 measures the size of the ellipse. Let, r 
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Then r is the sum of the semi-major and semi-minor axes. It is known 
[12, 16] that the convergence rate of the scheme is bounded by r“ N . Hence, 

as £ increases the approximation converges faster. £ is determined by 
the closest singularity. Suppose that this singularity occurs at x, "y, 
then 



II 

CM 

<** 

2(x 2 y 2 + 1 

♦/ 

— ? —9 9 —9 

(x + y -1) +4y ) 

. 


For 

f(x) = tanh(M(x 

- Xq)) we 

have 

that x — Xq and y 

tri 
2M * 

Thus, as 

x 0 

varies, y is fixed while 

X 

changes. It is easily shown 

that 

d 

a5 2 

(l 2 ) > 0. Thus, 

for fixed 

y> 

2 

£ is a minimum at 

o 

II 

lx 

and 


£ increases as |x| = |x^| increases. Hence, as xq approaches the 
boundaries, ±1, the rate of convergence increases. Also, as M increases, 
i.e., the function has a larger gradient, then y decreases and £ 
decreases and so the rate of convergence decreases. For interpolation 
approximations, both at Chebyshev and uniformly spaced points, a similar 
phenomenon occurs but the quantitative analysis is more complex [12]. 

For uniformly spaced collocation points in [0,1], the ellipses are 
replaced by the curves u(x,y) = constant where 

u(x,y) = 1 - xln/ x 2 + y^ - (1 - x)ln/(l - x)^ + y^ + y arctan — . 

x - x - y 

By examining graphs of this curve [11, p. 249] one sees that having the first 
singularity at x^ + i^ increases the size of the region as x^ moves 

toward the boundaries. As before, this increases the rate of convergence. 
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In Figure 6, we consider the same case as Figure 5a for the function 
f(x) = tanh[Q(x - xq)]. Here, Q is a function of M and Xq, specifically 


Q = 


xM 

2 


M 2 + 1 

1 + M 2 (1 - Xq ) 


M - 2, 4, 8, 15, 32. 


At the center, x Q = 0, Q = — while near the edge 

M 2 " 2 

Xq 21 1 and Q 21 ~~ 2 ~ • With this scaling the L z error is essentially 
independent of Xq* This indicates that an adaptive collocation method could 
be useful [2, 8]. 

In order to further investigate the resolving power of the schemes, we 
repeat the experiment of Figure 5 but for a function that is not analytic. In 
this case, our previous analysis is no longer valid. We choose 


f (x) 


! sgn(n) |n | > l 

| (3n 5 - lOn 3 + l5n) |n| < l 


(4.1) 


where n = M(x - x ). Hence, f(x) = -1 when x < - 77 , f(x) = +1 

(J U M 

when x > x~ + 77 and is a quintic polynomial in between. Furthermore, f (x) 
0 M 

has two continuous derivatives, but the third derivative is dicontinuous at 

x = * 77 • Thus, as before, x n denotes the center of the "jump” and the 

0 M u 

9 

gradient becomes larger as M increases. In Figure 7a, we plot the L“ 
error for Chebyshev collocation with 31 nodes. As Xq goes toward the 

boundary, there is a small decrease in the error, but not as pronounced in 
Figure 5a. Even more surprising is the fact that the decrease in error is 
greater for M = 32 than for M = 2. Thus, in contrast to Figure 5a, there 
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is no longer a sharp decrease for smoother functions as xq approaches 1 * 
When using uniformly spaced points the absolute error is larger than when 
using Chebyshev points. However, now there is a large decrease in the error 
as Xq approaches the boundary. We compare the case M = 16; for Chebyshev 
collocation the error decrease by about 2 orders of magnitude as Xq varies 
from the center to the edge. For uniformly spaced points the error decreases 
by about 6 orders of magnitude. This is despite the fact that the Chebyshev 
collocation method bunches the points near the edge. We also note that 
nothing special happens when Xq is sufficiently close to the boundary that 

the discontinuous third derivative at x. + 77 is no longer in the domain. 

U M 

In Figure 8 a, we study a similar phenomena. In this case, we study the 

2 

L error as we vary the strength of the singularity. We consider the 
function u(x) = H(x - Xq) * (x - Xq)^, where H(x) is the Heaviside 

function. Thus u(x) has a discontinuous M-th derivative. As expected, 
based on previous cases, we see that when the high order derivatives are 
discontinuous that the Chebyshev collocation method resolves the functions 
best when the discontinuity is near the boundary. However, when low order 
derivatives are discontinuous than the differential between the edge and the 
center decreases. For a step function, M = 0, the error oscillates with equal 
amplitude throughout the domain. As x approaches the boundary only the 
frequency of the oscillation changes. In Figure 8 b, we again see that the 
same qualitative picture occurs when the collocation is based on uniformly 
spaced points. We also see that global collocation based on uniformly spaced 
points is not convergent when the function is not smooth. This divergence is 
amplified if the discontinuity occurs near the center of the domain. In this 
case, the divergence is no longer caused by roundoff error. Rather it already 
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occurs at moderate values of N and begins at larger error levels. For f(x) 
= |x| it can be proved [16] that collocation based on uniformly spaced nodes 
converges only for the points x = 0, +1, -1. 

In order to further study the resolving power of the global schemes near 

x < 1 

x 2 1 

We plot the pointwise error in Figure 9a for both Chebyshev nodes and for 

uniformly spaced notes. For uniformly spaced nodes, the error is very small 

in the interior, (see Figure 9b for a logarithmically scaled plot) but is very 

large near, i.e., within 0(“-) , x = 1. From Figure 9b we see that the 

error is larger near x = -1 than in the center. For the Chebyshev nodes, 

the error is more global, but the large error near the boundary is confined to 

an interval of size 0 (”r) . 

N 2 

We next consider the partial differential equation 


the boundary, we consider the function 


f( 


■■ ■ 


u = u -1 <x< 1, t>0 

t x — — 


u(x,0) = f ( x ) u( 1 , t ) = g(t). 


We first discretize (4.2) in space using 


v = 
t 


(4.2) 


D v 


(4.3) 
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where D is the matrix derivative based on the collocation points xq,...xn 
and v is the vector of the dependent function evaluated at the collocation 
nodes. We shall further assume that the point x = 1 is a collocation 

point. As before D is explicitly given by (2.12). To advance (4.3) in time 
we use a four-stage fourth order Runge-Kutta formula. 

In studying (4.2) we shall be interested in both accuracy and stability 
properties of the algorithm. For stability we need to distinguish between 
space stability and time stability [6]. By space stability, we mean the 
behavior of the approximation v as the number of modes N increases when 
0 <_ t <_ T. By time stability we mean the behavior of v as time increases, 
for fixed N. Since, D can be diagonalized the scheme is time stable 
whenever all the eigenvalues of At-D lie in the stability region of the 
Runge-Kutta scheme. This does not necessarily prove space stability since the 
norm of the matrix that diagonalizes D depends itself on N. Obviously, the 
spectral radius of D and also the maximum allowable time step depends on the 
implementation of the boundary conditions. 

Since the temporal accuracy is lower than the spatial accuracy the 
maximum At allowed by stability considerations will not yield very 
accurate approximations. However, by decreasing the time step we can increase 
the accuracy of the solution. This general technique works equally well for 
nonlinear problems. When the model equation (4.1) is replaced by a more 
realistic system with several wave speeds then the stability limit will also 
give approximations that are accurate [1]. Also, when one is only interested 
in the steady state then frequently the time step can be chosen by stability 
considerations alone. An alternative, which will not be persued in this 
study, is to use spectral methods also in the time domain, e.g., [4, 21]. 
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In order to measure the accuracy of the approximation, we shall choose 
u(x,t) = f(x - t) for some f(x). Hence, the approximation can be compared 
pointwise with the analytic solution. The boundary data is then given by 
g(t) = f ( 1 - t). We shall measure the error either pointwise or else in a 
weighted Lr norm given by (3.2). 

We first study the effect of the boundary treatment on the stability and 
accuracy of (4.3). One property of global methods is that the approximation 
is automatically updated at all collocation points including the boundaries. 
Thus, if one wished, the scheme could be advanced without ever imposing the 
given boundary data; but this would be an unstable scheme. For a multistage 
time scheme, one can impose the boundary conditions at any stage one wishes. 
We now consider (4.1) with f(x) = sin(Trx). In Figure 10a, we impose the 
given boundary condition after each stage while in Figure 10b we impose the 
boundary condition only after the fourth stage. We define the Courant number, 
CFL, by 

CFL = N 2 At. 

In both plots, 10a and 10b, we display the error for several values of the 
Courant number. We see that imposing boundary conditions after each stage 
allows a larger maximum stable CFL number. For the four stage scheme, the 
maximum CFL is about 35. However, for smaller time steps the error is 
slightly larger than when one imposes the boundary condition only at the end 
of all the stages. One also sees that for a given error level that the 
approximate solution is essentially independent of the time step below some 
critical time step. As one demands more accuracy the necessary CFL number 


decreases. 


For a smooth solution, the necessary time step depends 
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exponent i ally on N. The largest stable At do not give accurate solution 
at any error level. We also found that the error grows in time if the 
solution is not sufficiently resolved in either space or time. There was no 
growth when N was large enough and At was sufficiently small. 

In Figures 11a and lib, we again plot the L 2 error for approximating 
(4.2) with f(x) = sin(x) as N increases and for a selection of CFL 

numbers. In this plot, we choose a different sequence of collocation points 
given by 


Xj = -(1 - a) cos gl + a(-l + ^-) j = 0, . . . ,N 


(4.4) 


so Xq = -1 and x N = 1. These points are a linear combination of Chebyshev 
nodes and uniformly spaced nodes. Letting 


a = 


(3 ~ 1)(1 - co&jj) 

2 , . IT -v 

ii - a - “V 


Then 


X (3 - 1) 1 


4N 


1 - 


4N 


(4.5) 


0(i) when 3 = 0(1), 


and we find that 


new spacing at edge 

Chebyshev spacing at edge 


(4.6) 
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We solve (4.1) by using the derivative matrix (2.12). We do not use a mapping 
to Chebyshev collocation nodes. In Figure 11a, we choose 6=2, i.e., a 

spacing at the edge twice as coarse as the usual Chebyshev spacing. We see 
that in this case we cannot increase the allowable time step beyond the 

stability condition for the Chebyshev nodes. Hence, the stability condition 
is not directly related to the minimum spacing. In Figure 11b, we display the 
error for 6 = V 2 > i.e., a spacing twice as small as the Chebyshev spacing 

near the wall. In this case the largest stable time step is reduced compared 
with the Chebyshev nodes. In this example, we have considered constant 

coefficients. For a problem with variable coefficients it is possible that 
coarsening the mesh near the boundary will allow a larger time step. This is 
because the coarser mesh near the boundary may just counteract the behavior of 
the variable coefficients near the boundary. 

In Figure 12, we consider uniformly spaced nodes, i.e., a = 0. From 

Figure 12, we see that even for small CFL levels that the error first 
decreases but then increases as N gets larger. These calculations were 
carried out in double precision on the CRAY. Nevertheless, it is difficult to 
distinguish between a mathematical instability and an instability caused by 
rounding errors on the computer. 

In Figure 13, we consider the differential equation 

u = -xu , -1 < x < 1 (4.7) 

t x — — 

u(x,0) = f(x). 

For this differential equation, we do not specify boundary conditions at 
either end of the domain. The solution is given by u(x,t) = f(xe t ) and 
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so u(x,t) decays to a constant value. It is easy to verify that the 

eigenfunctions of the spatial operator are v^(x) “ with corresponding 

eigenvalue A . = - j , j = 0,...,N (see also [5]). Hence, the stability 

C 

condition is At <_ — where C depends on the details of the explicit time 
integration scheme. Since u(x,t) is almost constant for large time the 

error levels become very small. In figure 13, we plot the error at 

time t * 10, with f(x) ■ sIotx. For the fourth order Runge-Kutta scheme 

C ~ 2.8 is the stability limit. In Figure 13a, we use double precision on 
the cray (about 30 significant figures) while in Figure 13b we only use five 
significant figures. We define the CFL number for this problem as 

CFL = NAt. 

Comparing 13a with 13b, we note that CFL = 2 is stable using double 

precision but is unstable using only five significant digits. The effect of 

roundoff on stability is studied in [24]. For this case the effects of 

roundoff are important only for time steps very close to the stability 

limit. The effect of roundoff is more pronounced when At ~ 1/N than 
2 

when At ~ 1/N . 

We further see from this case that the maximum allowable time step is not 
necessarily related to the minimum spacing in the grid. In this case, the 
fact that no boundary conditions were specified allowed At to vary with 
1/N rather than the usual 1/N . We also saw a similar phenomenon where 
coarsening the mesh near the boundary did not allow a larger maximum time 
step. A similar conclusion was found by Tal-Ezer [22] for the Legendre-Tau 
method which has a time step limitation that depends on 1/N even though the 
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minimum grid spacing is 1/N . Thus, to find the stability limit, one must 
analyze the derivative matrix appropriate for each case rather than using a 
heuristic approach based on the spacing between collocation nodes. 

We finally discuss the solution to the Poisson equation (3.5). The right 
hand side and boundary conditions are chosen by deciding a priori on the exact 
solution. We solve (3.5) by a Chebyshev collocation method in each 
direction. The matrix equation that results is solved by a multigrid 
technique [25]. 

In Figure 14 we consider the case where the exact solution is u(x,y) = 
simry tanh(M(x - x^)). Thus, the solution is smooth in y and has a 
gradient in the x direction. The sharpness of the gradient and its location 
are given by M and Xq respectively. Hence, this models boundary layer 
type behavior. As before (see Figure 5a) when M is not too small then the 
approximation is more accurate when the gradient occurs near the boundary. 
For sharp gradients, i.e., M very large, the gradient is not resolved by the 
mesh and the Chebyshev L error is approximately independent of the position 
of the gradient. As shown by Figure 5b, this increased accuracy in the 

boundary layer is not only due to the increased number of collocation points 
in the "boundary layer". Rather it is due to properties of global 
approximation techniques. It is of interest to note that for M = 1024, i.e., 
a discontinuity, that the error is almost constant. However, for M = 64 and 
256, i.e., a sharp gradient, there are peaks in the error as xq approaches a 


collocation node. 
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5. CONCLUSIONS 

We consider the properties of global collocation methods to problems in 
approximation theory and also partial differential equations. In particular, 
we study concepts that have been used by many authors without verification. 

In order to be able to study differential equations for a general 
sequence of collocation nodes we calculate the approximate derivative by a 
matrix times vector multiply. For Fourier or Chebyshev methods one could also 
use a FFT [5]. For a small number of nodes, N ~ 64, the matrix multiply is 
faster than the FFT. For sufficiently large N the FFT is always faster 
since it grows as NlogN rather than N . The exact cross-over point between 
the two techniques is very machine dependent as well as software dependent. 
There obviously are differences between scalar, vector, and parallel 
computers. Nevertheless, for practical N used in most partial differential 
equation solvers the matrix multiply is not much slower than the FFT. Hence, 
we only use the matrix multiply technique due to its greater generality and 
flexibility. 

It follows from the results presented in Section 4 that a global 
collocation method must be distinguished from a local finite difference or 
finite element approximation. In particular, the greater density of points, 
for a Chebyshev collocation method, near the boundary does not give increased 
accuracy, for a smooth function, near the boundary. The extra density near 
the boundary is needed to counteract the tendency of polynomials to have large 
errors near the edges of the domain. 

Chebyshev collocation methods do have lower errors when sharp gradients 
are near the boundary than when they are in the center of the domain. Similar 


results occur when there is a discontinuity in some derivative. However, 


-28- 


qualitatively similar results are obtained using uniformly spaced nodes. 
Thus, the increased resolution near the boundary is due to the global nature 
of the approximation and not the bunching of collocation nodes. Of course, in 
terms of absolute error, it is preferable to use Chebyshev collocation rather 
than uniform collocation. This indicates that domain decomposition methods 
should be advantageous [9, 10] but not for shocks. In fact, even in cases 

where collocation based on a uniform mesh should converge the actual 
interpolation process on a computer eventually diverges due to roundoff 
errors. These roundoff errors contaminate the results for relatively small 
N. 

As a further distinction between global and local techniques we consider 
the aliasing limit. For a Fourier (periodic) method we need 2 points per wave 
length to resolve a sine wave. For a Chebyshev method we need tt points 
per wavelength. The difference between 2 and tt is not due to the 
different distribution of points in these techniques. Polynomial collocation 
based on uniformly spaced points again needs tt points per wavelength. 
Furthermore, for other functions, e.g., tanh x, one does not observe any sharp 
aliasing limit. Thus, one can not speak of number of points per wave length 
for general functions on nonperiodic domains. 

An alternative to improving the accuracy of an approximation is to map 
the x domain [-1,1] onto another computational domain s, for simplicity 

again [-1,1]. The above conclusions do not extend to such mappings. First, a 
polynomial in s is no longer a polynomial in x. Hence, in the physical 

space x we are not considering polynomial collocation methods. In addition, 

2 2 
the L norm in s-space corresponds to a weighted L norm m x-space. 

Hence, it is difficult to measure the effectiveness of such mappings. In 

practice [2] has shown that in some cases adaptive mesh mappings can be 

effective for spectral methods. 
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The results obtained for approximating solutions to elliptic partial 
differential equations seem to correspond to the results for the approximation 
problem. Again one can not interpret the properties of a Chebyshev 
collocation method in terms of finite difference properties. Such concepts as 
number of points in a local region are not meaningful. If one chooses another 
set of collocation points then there are two ways of implementing the 
method. One can map one set of points to the other and then use a Chebyshev 
method in the computational space. This introduces metrics into the 
equation. Alternatively, one can solve the equation in physical space using 
the general derivative matrix (2.12). We have not investigated the 
differences between these two approaches. 

For a time dependent partial differential equation, the study is more 
complicated. First, there is an accumulation of errors as time progresses. 
Thus, for example, for a stationary problem one can distinguish between the 
discontinuity being at a node or in between nodes. For a time dependent 
problem the discontinuity is moving and so all effects are combined. This is 
especially true for systems with variable coefficients where there is coupling 
between all the components. 

Also, there is the question of stability in addition to accuracy. Thus, 
we have found that the implementation of boundary conditions influences both 
the maximum time step allowed and the accuracy. At times an implementation 
which increases the stability will decrease the accuracy. 

We also found that there is no direct correlation between the smallest 
distance in the mesh and the maximum allowable time step. Coarsening the mesh 
near the boundary does not allow a larger time step. This again demonstrates 
the fallacy of describing a global method in terms of local behavior. As is 
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well known, for wave equation type problems one should not choose the maximum 
allowable time step allowed by stability. Since we use a fourth order 
accurate method in time but a spectrally accurate method in space one should 
choose a smaller time step to compensate. Thus to achieve time accuracy there 
is no need to increase the 0(1/N^) time step restriction for hyperbolic 
equations. For stiff problems or if one is not interested in time accuracy 
then one may wish to exceed the stability restriction. Furthermore, for 
parabolic equations At 0(1/N^) which is much too restrictive. As 
before, one can consider other sets of collocation points. Again using 
mappings or the derivative matrix based on these nodes give rise to different 
schemes. 
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APPENDIX A 


In Section 2, we saw that given the collocation points xq,...,x n and 

N + 1 functions Uj(x) then the derivative matrix, D, is determined by 

demanding that D times (<}>j (x Q ) , . . . ,<|k (x n > ) C give the exact derivative at 

the collocation points. Let D = (dj^) and define the matrix U by = 

Uj(x^) j,k = 0,...,N. Given the matrices D and U we denote the j-th 

column of these matrices as dj and Uj . Then each column of D is 

determined by the equation 

Ud. = u : 

J 3 

(Al) 

at all collocation points x^, k = 0,...,N. 

If we wish D to be exact for M > N + 1 functions, then in general there is 
no solution. Instead we can demand that D give the smallest L error over 
these M functions. Intuitively if D is almost exact for many functions, 
then it should be a good approximation to the derivative. In particular, one 
may choose functions that are more appropriate to a given problem than 
polynomials. Choosing D to give the least sqares minimization is equivalent 
to demanding that 


u fc ud. = u T u: 

J 3 


(A2) 


instead of (Ai). It is easily to verify that 
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and 


(D u) jk ■ J 0 u i (x j )u i (x k > 


(U t u') = l u (x. )u'(x ), j,k = 0,...,N. (A3) 

J ^ i=0 1 K 1 J 


We now define 


M 

v k (x) = l u i (x k )u i (x). 


k = 0, . . . ,N 


(A4) 


then 


v'(x.) ■ ‘"‘“JV 

Assuming det U ^ 0 then the v^Ct) are linearly independent. It also 
follows that D is exact for the N + 1 functions v^Cx) at the collocation 
points. Hence, demanding least square minimization for u^Cx), i = 
is equivalent to demanding exactness for v^x), i = 0,...,N given in (A4). 

We next extend this by letting M become infinite and replacing the sums 
by integrals. Thus, given the continuum of function u^(x) and demanding 
that D be the best least squares approximation to the derivative at the 
collocation points is equivalent to demanding that D be exact for the 
N + 1 functions 


m 

v k (x) = / u i (x k )u i (x)di. 


(A5) 
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To demonstrate this, we consider a specific example. Let u^Cx) be the 
functions sin(kx) and cos(kx) for 0 < k < and choose N + 1 
uniformly spaced collocation points, x^ . Since k £ -jp- we are always below 
the aliasing limit. It follows from (A5) that 


Nil 

2 


v„(x,) = / 
0 


'k' J 



. Nir , . 

sin — (x^ - x^) 

x. - x 
J k 


Nir 

2 


j * k 


j = k 


(A6) 


These functions, Vj(x) are ^ nown as SINC functions and have been used for 
interpolation formulae [19]. Demanding that D be exact for Vj(x), j = 
yields the derivative matrix 


d 


jk 



J * k 


j - k 


which is an antisymmetric matrix. We also note that this matrix resembles the 
derivative matrix for the Chebyshev nodes [7]. 
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APPENDIX B 


In this section, we present the proof that Chebyshev collocation at the 
standard Gauss-Lobatto points is stable for solving scalar hyperbolic 
equations. This result was given in [7] without proof. 

Consider the collocation points 

Xj = cos , j = 0,...,N. (Bl) 

Let u be the solution to 


u t = u x u(l,t) = 0 

u(x,0) = f(x). 


(B2 ) 


If v is a N-th order polynomial which is found by collocation at Xj (see 
[5], [7]), then v exactly solves the modified equation 


v 


t 


v + 
x 


o + x)t nV t) 

N 


v(l,t) = 0 


(B3) 


where 


N 


d 

dt a N 


and 


N 

v(x,t) = l a k (t)T k (x). 


(B4) 


We need the following fact: 
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Lenma: 


U2a N " a N-l )2] = " 4N(2a N - VW* 


(B5) 


Proof : 


^ 2a N " a N-l^ ^ 2 ^ 2a N a N-l ^ 2 


da N __ da N— 1 > 
dt “ dt ' 


(B6) 


comparing the coefficient of i n (S3) we find that 


da N-l da N 

= 2Na„ + 2 


dt 


N dt 


or 


da da 

2 _JL - — Sti = -2Na M . 
dt dt N 


Inserting this into (B6) gives the lemma. With this lemma we prove the 
following theorem. Let 


1 


j = 0,N 


Cj (2 j * 0,N 


then 


Theore®: Let v solve (A2.3), then if 


3FT< 6 iT iisa 
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d^_ 

dt 


{|n l T: (1 + X j Kl “ 6x.)v 2 (x.,t) +|g- <2^ - a^) 2 } < 0 (B7) 


J J 


and so the solution v is stable. 


Proof: We multiply (B3) by (1 + x^)(l - gx,. )v(x. ,t) and sum over 

NUj J J J 

the collocation nodes. We then have 


a If | ^ (1 + x j ,(1 ‘ Sx j )v j ' ¥ | it 1(1 + V (1 ‘ 6x j^ vv x 


+ + V 2<1 -eVVV v< V' 

J J 


However, the last term is zero since T' (x^ ) =0 at interior points, 
1 + Xj = 0 at Xj = -1 and v(xj) =0 at Xj = +1. Furthermore, if 


(B8) 


f(x) = l bj Tj (x) 


f£P 


2N-3 


then 


I ^ 77 ' 

3 


f(x) 


/ 1 - x 2 


dx + xb 


2N* 


By 

algebra, it can 

be 

verified 




(i 

(1 

+ x) (1 - 3x)vv x 

is 

b 2N = 

a j 

are the Chebyshev 

coefficients 


(i - e)Na; 


N 0 x 2N-K 
- " e(~ 7 ! — )a XT a v 


4 N N-l 


where as before 
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kklh" + V (1 


1 (1 + x)(l - 0 x)vv 


J J 


SX.W. ■/, 


dx 


/ 1 - 


(B9) 


+ 1 id - - t^vn 1 


Integrating by parts and using the fact that v(l,t) = 0 we find that 


TT d_ y 1 

2N dt : c, 
J J 


(1 + X . ) ( 1 - 6X >V 2 (X ) = “y / 


1 (1 - 0 - 0 x + 0 x 2 )v 2 (x,t)dx 


-1 


(1 - x) 


/ 1 - 


(BIO) 


+ 2 ^ ®^ Na N ^ 2 


Using the lemma this is equivalent to 


Se fa j ^ u + V (1 ' 6x 3 )v2<x 3’ t> + W <2a » ' Vi>*» 


(Bll) 


- I f 1 3 ~ gx ±J±- v 2 (x,t)dx - 1 [(30 - 1)N - 0]a 2 . 

- 1 - - X) /7T7 


(1 


If B £ 3 - , then the integral term is negative while if 0 >. ~ -j , 

the second term on the right hand side is also negative. Hence, when 


then 
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— < R < — then the right hand side of ( B 1 1 ) is negative and the theorem 
3N-1 — p - 5 


is proven. 


If we choose the special case 8 = — then (Bll) becomes 


Ir & 2 t : (1 + x j )(1 -f V^v 0 + w (2a N ■ 3 n-i )2 ^ 


(B12) 


" "To/ 


10 i . (1 - ~ v 2 (x,t)dx - |q( 7N - 4)a 2 . 


1 ( 1 - x)/l - X 2 


As a corollary, this theorem implies that all the eigenvalues of D lie in 
the left half of the complex plane. 
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Figure lb. Pseudospectral Chebyshev approximation to 


sin(20x - M), 0 < M < 2- , 


with N - 29 
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Figure 3a. Pseudospectral Chebyshev approximation to sin(MTrx). The 

different graphs represent M = 2, 4, 8, 16, and 32. The 

error is plotted as a function of — with 

M 1 ' 


several N 
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Figure 3b. Pseudospectral Chebyshev approximation to tanh(Mx) for M = 
2, 8, 16, 32. We plot normalized L error as a function 


of N/M for several M. 
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Figure 4. 



12 ERROR IN 


-49- 



Figure 5a. Pseudospectral collocation with 31 nodes for f(x) = 

tanh(M(x — Xq)) with M = 8, 32, 512, 2048. xq varies 
between the center, x Q = 0, and the edge, x Q = 1 . 
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Figure 5b. Same case as Figure 5a but using uniformly spaced collocation 
points. Now M = 2, 4, 8, 16, and 32. The error is the 

same Chebyshev norm as in Figure 5a. 
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Figure 9b. The uniformly spaced nodes of Figure 9a plotted on a 


logarithmic scale 
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Figure 10a. Pseudospect ra I approximation to (4.1) with f(x) = sinirx. 

A four-stage fourth-order Runge-Kutta formula is used and 
boundary conditions are imposed after every stage. Each 
graph represents a different time step, i.e., CFL number with 

n 

an increase of /2 between graphs. The Lr error at t = 
1 is given as a function of N. 





Figure 10b. Pseudospectral approximation to (4.1) with 

f(x) = sinxx. A four-stage fourth-order Runge-Kutta 
formula is used but with the boundary condition imposed only 
once after the completion of the four Runge-Kutta stages. 
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Figure 11a. Pseudospectral Chebyshev approximation to (4.1) with 

f(x) = simrx, 4 applications of the boundary conditions 
and 3 = 2, i.e., mesh is twice as coarse as a Chebyshev 
grid near the boundary. Each graph represents a different 
CFL number, increasing by a factor of /2. 
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Figure lib. Pseudospectral Chebyshev approximation to (4.1) with 

f(x) = sinrrx, 4 applications of the boundary conditons 
and g = , i.e., twice the density of Chebyshev spacing 
near the boundary. Each graph represents a different CFL 
number, increasing by a factor of /2. 
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Figure 12. Same as Figure 10a but with uniformly spaced nodes. Based 
on double precision on the CRAY-XMP. 
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Figure 13b, Same as Figure 13a but using only 5 significant digits. 
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Figure 14. Pseudospectral Chebyshev approximation to the solution of a 

Poisson equation. The exact solution is u(x,y) = sinny 

tanh(m(x - Xq)) with M = 4, 16, 64, 256, 1024 and N = 

o 

17 modes in each direction. We plot the L Chebyshev 

error as a function of x n . 
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